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Development of a general aerodynamic model that is adequate for predicting the forces 
and moments in the nonlinear and unsteady portions of the flight envelope has not been 
accomplished to a satisfactory degree. Predicting aerodynamic response during arbitrary 
motion of an aircraft over the complete flight envelope requires further development of the 
mathematical model and the associated methods for ground-based testing in order to allow 
identification of the model. In this study, a general nonlinear unsteady aerodynamic model is 
presented, followed by a summary of a linear modeling methodology that includes test and 
identification methods, and then a progressive series of steps suggesting a roadmap to 
develop a general nonlinear methodology that defines modeling, testing, and identification 
methods. Initial steps of the general methodology were applied to static and oscillatory test 
data to identify rolling-moment coefficient. Static measurements uncovered complicated 
dependencies of the aerodynamic coefficient on angle of attack and sideslip in the stall region 
making it difficult to find a simple analytical expression for the measurement data. In order 
to assess the effect of sideslip on the damping and unsteady terms, oscillatory tests in roll 
were conducted at different values of an initial offset in sideslip. Candidate runs for analyses 
were selected where higher order harmonics were required for the model and where in- 
phase and out-of-phase components varied with frequency. From these results it was found 
that only data in the angle-of-attack range of 35° to 37.5° met these requirements. From the 
limited results it was observed that the identified models lit the data well and both the 
damping-in-roll and the unsteady term gain are decreasing with increasing sideslip and 
motion amplitude. Limited similarity between parameter values in the nonlinear model and 
the linear model suggest that identiliability of parameters in both terms may be a problem. 
However, the proposed methodology can still be used with careful experiment design and 
carefully selected values of angle of attack, sideslip, amplitude, and frequency of the 
oscillatory data. 
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Multiple Correlation Coefficient 
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non-dimensional time constant, r, = V l(!'b [ ) ; longitudinal case i = 1, lateral case i = 2. 


v = measurement noise 

co = angular frequency, rad/sec 
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A = amplitude (subscript) 

E = measurement 

Aerodynamic Derivatives 

( 00 ) = C a ^ = where a = (C D ,C L ,C Y ) or ( C m , C t , C n ) and £ = (a, 7 ?) 


C %<“)- C -fe 


= 0Q 


where a = (C D , C L , C r ) or (C m , C ( , C n ) and | 2 = (p, q, r) 


I. Introduction 

D EVELOPMENT of a general aerodynamic model that is adequate for predicting forces and moments in 
nonlinear and unsteady portions of the flight envelope has not been accomplished to a satisfactory degree. 
Aerodynamic phenomena such as shock waves, separated flows, and vortex flows can cause faulty prediction of 
aircraft response when conventional linear aerodynamic models are used. Predicting aerodynamic response during 
arbitrary motion of an aircraft over the complete flight envelope requires further development of the mathematical 
model and associated methods for testing that will allow identification of the model. 

Development of the aircraft aerodynamic model began in the early days of flight. Byran’s formulation 1 in 1911 
established the conventional assumption that aerodynamic forces and moments depend only on instantaneous values 
of aircraft states and can be expressed in a linear polynomial expansion. These assumptions work well, depending on 
the particular configuration, but are primarily restricted to linear portions of the flight envelope. Adding nonlinear 
terms to the series expansion 2 extends the capability of the model allowing a certain class of model nonlinearity to 
enter the response but still no memoiy of previous states is captured. Improving this formulation to accommodate 
unsteady responses requires additional differential equations to capture the lag in response. Research efforts 
addressing the nonlinear unsteady problem began in the 1970’s with recognition that the Volterra functional was a 
desirable mathematical structure 3 . This model structure, in the form of indicial models, is sufficiently general to 
capture a large class of nonlinear unsteady behaviors. A number of researchers have investigated this modeling 
problem. The modeling of unsteady aerodynamic characteristics of an aircraft from experimental data was first 
reported in Ref. 4. The model equations were formulated as state-space equations. Further improvement to this 
modeling and parameter estimation followed in Refs. 5-7. Many other researchers have approached the modeling 
problem using an indicial function formulation. In Refs. 8-10, linear models were applied to wind tunnel data from 
different testing facilities for different aircraft models. In Refs. 11 and 12, nonlinear indicial response models were 
applied to the rolling 65-degree delta wing and in Ref. 13 to the prediction of a dynamically stalling wing. Nonlinear 
indicial models were also used in model identification from wind tunnel oscillatory data using least squares and 
maximum likelihood principles 14 . In Ref. 15 a methodology was introduced that allowed efficient identification of a 
linear unsteady model for an F-16XL configuration in planar forced oscillation. In this approach the aerodynamic 
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model was formulated with static and steady-flow damping terms separate from unsteady indicial components. 
Enhanced identification and efficient testing was accomplished by using wide-band inputs. In Ref. 16 additional 
experiments were performed to validate the methodology, mathematical model, and test techniques used in Ref. 15. 

The present study extends the work in Ref. 16 as part of the ongoing effort at NASA Langley Research Center 
(LaRC) to develop a more general aerodynamic model and associated test methods. To date, development and 
application of the indicial model methodology at LaRC has been done with a number of limiting assumptions. This 
is a natural result of the research progressing from specific simpler models to more general complex models. In this 
paper, a general nonlinear unsteady aerodynamic model is presented, followed by a summary of a linear modeling 
methodology that includes test and identification methods. Then a progressive series of steps that remove limiting 
assumptions and suggest a roadmap toward development of a general nonlinear methodology are presented. 
Included in this paper are results from application of two initial steps of the general methodology applied to static 
and oscillatory test data to identify rolling-moment coefficient. 

The general methodology for testing and data analysis is being developed using data from both wind tunnel and 
water tunnel tests. Results presented in this paper are taken from water tunnel experiments. This choice provides 
data at relatively low cost and does not affect conclusions drawn regarding the methodology. Tests are done, under a 
Phase III SBIR with NASA LaRC, using an advanced dynamic test rig at Rolling Hills Research Company. 

II. General Aerodynamic Model 

Assumptions used for this general model impose only a few limitations. For this case only rigid-body responses 
are considered and relatively constant flight path properties (e.g., density, viscosity, dynamic pressure) are assumed. 
And assuming that the aerodynamic coefficients are single-valued functions of aircraft state variables a, (1. p, q, and 
r, then according to Ref. 3, these coefficients can be formulated as 


C a (0 = C fl (0) + } C (t - r; £(r )) r -f^dr 
o bl dr 

k o a dr 


(1) 


where C a (t) is a coefficient of aerodynamic force or moment, CJO) is the value of the coefficient at initial steady- 
state conditions, and C a d(t) is a vector of indicial functions whose elements are the responses in C„ to unit steps in q. 
E, is defined as 


I = [£i :£>] = [« P'-P 1 r f 

and i is the characteristic length ( i = c / 2 or i = b w / 2 ). The indicial functions approach steady-state values with 
increasing values of the argument ( t-r ). To indicate this property, each indicial function can be expressed as 


c ap .(t-T’ 4dz-)) = C . (oo; £(r)) - F ( t - r;4Tr)) 

f/ f/ 


(2) 


where C ac (oo ;^(r)) is the rate of change of the coefficient C a with A- , while the remaining variables c arc fixed at 
the instantaneous values of c(t) , and F is the deficiency function. This function approaches zero for 

f/' 

(t-r)— > oo. When Eq. (2) is substituted into Eq. (1) the terms involving the steady-state parameters can be 
integrated and Eq. (1) becomes 
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( 3 ) 


C a (0 - C a (co; m -] F a At- r; £(r )) r £ (r)rfr 
o dr 

~~\F (t-r,g(T)f -^-^ 2 (T)dr 

V o 2 dr 

where C a (oo; c(t)) is the total aerodynamic coefficient that would correspond to steady flow with c fixed at the 
instantaneous values c(t). F a and F a ^ are vectors of deficiency functions. Further simplification of Eq. (3) is 

achieved by assuming that the coefficients C a are linearly dependent on the angular rates c 2 • Applying a Taylor 
series expansion about = 0 to terms in Eq. (3) and keeping only the linear terms results in 


/ fir 1 

C a (t ) = C a (o°;£ 1 (0,0) + -— H«>-,m,0) r 52 

v d%2 

-\F (t-rMr), Of ^(T)dT 
o 1,1 dr 


(4a) 


or in simplified notation 


c a (0 = C a (oo; a (t),fl(t)) + y C ap (oo -a, P)p{t ) 
+ y C a q (°° -,cc,P)q+^C ar (OO \a,P)r(t) 

t 

~\F (t-T\a{T),p{r))d{T)dr 

o 

-J F an (t - T-a(T),P(r))P(T)dr 
o F 


(4b) 


where 


C„ 


_ dC c 




= and C a 

e <FL 


2V 


SC a 

~ rb 
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The deficiency function in Eq. (4) is postulated as 


F(t;a,P) = h(t;a,P)a(a,P ) (6) 

h(t;a,P) = e~ b(a ’ /})t (7) 

where a(a,P) and b(a, ft) are polynomials in (a, ft). When the differential operator, D = d/dt, is introduced, then the 
first convolution integral in Eq. (4) can be expressed as 
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(8) 


I F a a (* ~ *"> a(T),P(T))d(r)dT = D a(t) = rj a (. t ) 

o D + b { {a,fi ) 


and similarly for the second convolution integral 


D + b 2 (a,fi) 


DP(t) = Vp(t) . 


(9) 


Equations (8) and (9) represent first-order linear differential equations 


V a P) + b \ (a,/3)ri a (t) = a x (a,P)a 


(10) 


and similarly for rjp(t) . Substituting Eq. (8) and Eq. (9) into Eq. (4) results in 


C a (t) = C a (oo -a, P) + — C ap (go; a, P)p(t) 

+ y C Uq (oo; a,P)q(t) + y C ay (oo; a,P)r(t) . 
° l<a - /!| -Da(,)- Dfi V ) 


D + b { (a, P) 


D + b 2 (a,p ) 


(11) 


TIT. Linear Model Identification Methodology 

Linear model equations can be obtained from Eq. (11) by considering small deviations in state variables from 
steady-state trim conditions. Then C L ,C D , and C„, can be expressed as linear functions of (a.q) and 
C Y ,C f , and C„ as linear functions of ( p. p,r) . The resulting equation for the longitudinal aerodynamic coefficient 
is 


C a =C a {0) + C aa a + yC aq q 


-\F aa (t-T)CC{T)dT 

0 


and for the lateral coefficient 



~\ F a B (t-T)p(r)dr 

o p 


where 


(12a) 


(12b) 
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F a a 

F ap = « 2 e ”* 2 ' 

Rewriting Eq. (12a) with the state space form suggested by Eqs. (8) and (10) the linear state space form is 

ij(t) = + ad 

c a (0 = C aa (oo ) a {t) + y C aq (oo )q(t) — 77(0 ' 

The four unknown parameters, C aa (oo), C a ( oo), a, A , in Eq. (13) can be estimated, in principle, from 

measurement time histories of C a ,a,a and q. If these data are obtained from a one-degree-of- freedom oscillatory 
test in the wind tunnel (or water tunnel), then a = q and the aerodynamic model equations can be expressed in the 
form of a transfer function 


CJs) _ As 2 +Bs + C 
a(s ) s+b { 


(14) 


where s is the Laplace transform parameter and 


A 


~C a (oo) 

y a q 


B^C^W-a + ^yC^ (oo). 
C = b,C aa {K) 


(15) 


There have been several estimation techniques suggested for parameter estimation from measured data. In Ref. 
16 three different methods were used, 2-Step regression in the time domain, regression and maximum likelihood in 
the frequency domain. In the 2-Step regression parameter /?, is estimated first and then used as a known parameter in 
estimation of the three remaining parameters. The cost functions for the linear regression and maximum likelihood 
methods, respectively, have the forms 


and 


1 

J(0) = XI^O'XA +jco i ) + (Acof-C-jBco i )a(i) |" 
1=1 


(16) 




c a ( o- 


Acof - C - jB(o l 

A +M 


am 


(17) 


where N is the number of data points. Parameters are then obtained from Eq. (15). 

Oscillatory tests using a simple harmonic input are usually repeated at different frequencies. If the data are to be 
used to estimate (four) unsteady model parameters then six or more frequencies are recommended for better 
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statistical results. In order to avoid the large number of inns, the use of a wide-band input is recommended 15 . 
Specifically, the Schroeder sweep 17 was selected with specified amplitude to provide a flat power spectrum over a 
specified frequency range. Transforms of the time histories of the input and aerodynamic coefficients to the 
frequency domain were accomplished using a Discrete Fourier Transform (DFT) algorithm 1 6 . This algorithm utilizes 
a zoom transform and allows the transform to be performed over the frequency range corresponding to the wide- 
band input. 


IV. Nonlinear Model Identification Methodology 

A nonlinear modeling methodology that defines the mathematical model, algorithms for identification, and 
required test techniques is presented in this section only as a first candidate toward developing a more complete and 
general nonlinear methodology. Three steps in the candidate methodology suggest a progression for developing a 
general methodology. The modeling philosophy behind the progression is to utilize the least complex model that 
provides an adequate representation of the aerodynamic response; basically an application of the principle of 
parsimony. With this philosophy, a progression is made toward more complex models and correspondingly more 
complex experiment design. The simplifying assumptions associated with the least complex model are removed in 
progression, as required to achieve model adequacy. Before progressing through a series of dynamic tests, a 
complete static model is required to identity regions where static nonlinearities exist and to define the static term for 
the model. 

Step 1. In the first series of dynamic tests, a linear model methodology is applied. Dynamic tests begin with 
single-frequency, one-degree-of-freedom, small-amplitude motions including constant nonzero sideslip during 
longitudinal motions and nonzero mean sideslip during lateral-directional motions over a large range of angle of 
attack and sideslip. Inherent in these tests is the simplifying assumption of superposition, so longitudinal and lateral 
axes are treated independently. Using harmonic analysis 19 complexity of the aerodynamic response is assessed to 
determine where the least complex model can be applied with sufficient adequacy or equivalently to determine 
which simplifying assumptions are violated, if any. In general, it is expected that with small-amplitude maneuvers 
only a conventional linear response model or a linear unsteady model will be required and therefore only a 1 st order 
harmonic model (conventional in-phase and out-of-phase coefficients) will result from harmonic analysis. If a linear 
model is adequate then the linear model methodology applies and can be used to identify the full linear unsteady 
model from either single-frequency 10 or wide-band forced-oscillation tests 16 . If the linear model is not adequate then 
a more general model structure must be utilized. One indicator of model adequacy is the Multiple Correlation 
Coefficient (7?"). The metric value varies between 0 and 1 and is reduced from 1 as the output signal fails to be 
represented by the model. This metric can be computed as part of the parameter estimation process and can be 
expressed as 


l[z(i)-y(i)] 2 

R 2 = 1 — , 0 < 7? 2 < 1 ( 20 ) 

Y.[zU)-zf 


where 


and y(i) is the computed value of the measurement z(i). An important motivation for this step is that it reveals the 
functional relationship (for small amplitude motion) of the two dynamic terms with respect to angle of attack and 
sideslip. In general, it is expected that a nonlinear relationship exists for the dynamic terms as functions of angle of 
attack and sideslip. The small amplitude analysis captures the locally linearized model variation with angle of attack 
and sideslip. 

Step 2. In the second series of dynamic tests, a nonlinear model methodology may be required. These dynamic 
tests begin with single-frequency, one-degree-of-freedom, large-amplitude motions with constant nonzero sideslip 
during longitudinal motions and nonzero mean sideslip during lateral-directional motions over a large range of angle 
of attack and sideslip. With the assumption of small amplitude motion removed, the possibility of nonlinear 
response is more likely. As in the previous step, the simplifying assumption of superposition is still used so 
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longitudinal and lateral axes are treated independently. Although each axis, individually, may exhibit nonlinear 
response this assumption is equivalent to assuming that no aerodynamic coupling exists between axes. Again 
harmonic analysis is applied to determine the complexity of the aerodynamic response. The existence of higher 
order harmonics directly indicates the presence of nonlinear response in the model, although it does not indicate 
which terms are contributing to the nonlinear response. However, identification is possible since the static term is 
completely known from direct tunnel measurements and the unsteady term can be separately estimated using 
oscillatory coning experiments 16 ' 20 or through planar oscillation experiments 21 in heave and sway. The steady-flow 
damping term can be estimated knowing the other two terms or by using a curved- flow tunnel 21 ' 22 the term can be 
directly measured. 

A constraint for this study is that only conventional single-frequency, one-degree-of-freedom, forced-oscillation 
test data are considered. Under this constraint some analysis is still possible. If the static and steady-flow damping 
terms are assumed to be nonlinear and well represented by higher order polynomials, then the order determined by 
harmonic analysis also implies the order of the polynomial required. For example, consider a conventional harmonic 
model for rolling moment analysis 


AC, 


- C efjP + f C (f)P + f C (pP + 



CfpP 


(21) 


In this case, assuming (f> = <j) A sinatf and ft = cfcsvnd , it follows that a 1 st order harmonic model is 


2 

AC, =(j> A {C, sin 0 -k"C,p) sin cot + <f> A k(C, p + C,p sin 6) cosot 


(22) 


If a cubic nonlinearity exists in the static term, such that Eq. (21) contains the term C 3 ft , then this term produces 

a cubic-harmonic sine term. Using the identity, sin 3 cot = 1 / 4[3 sin cot - sin 3<ot 1 , it can be shown that both 1 st and 3 rd 
harmonic terms are added to the model to produce 


AC, = <]) A (C,p sin 0-k 2 C f p +(3/ ^)</ A C, 3 sin 3 #) sin 

+<j> A k(C fp +C,p sin#)cos®t . (23) 

-(1/4)^Q 3 sin 3 #sin3®t 


Eq. (23) shows that the cubic nonlinearity expresses itself as both a 3 rd order harmonic and an additional component 
added to the 1 st harmonic term. A similar analysis can be done for harmonics of different order and the same 
analysis can be applied to the steady-tlow damping term. Unfortunately, the unsteady term does not lend itself to 
this type of analysis. Thus, the presence of a 3 rd order term in the static data or in the steady-tlow damping term will 
result in a 3 rd harmonic during a general harmonic analysis, however, to determine the source of the higher harmonic 
response independent information about each of the three model terms is required. 

One possible approach to the identification problem is to assume that the steady-flow damping term is known 
approximately from Step 1. Given this approximation and the static measurements, identification of the unsteady 
term and the steady-tlow damping term can proceed using a relaxation method 14 . Since motions are limited (in this 
step) to one degree of freedom, forced-oscillations, a simplified version of model Eq. (11) can be used for model 
structure determination and parameter estimation. For examples in this paper, during large amplitude, body-axis roll 
oscillations ((/>a < 30°,) with a zero sideslip offset angle, the angle of attack will vary by less than 4 degrees. Even 
with a 15 degree mean offset in sideslip, the maximum angle of attack variation during oscillations is limited. 
Oscillations in a and /i during roll oscillations, vary according to the relations 

8 


American Institute of Aeronautics and Astronautics 



P = sin 1 {sin(fzi) sin(<9 0 ) cos(^ 0 ) - cos(^) sin(^/ 0 )} 


a = tan 


jcos(^) tan(<9 0 ) + 


sin(^)sin(y/ 0 ) 

COS^oJcOS^q) 


(24) 


where 9 0 and y/n are the nominal pitch and yaw offsets, respectively. In either example for this paper, it is reasonable 
to assume that C e can be expressed as C f {/3,p) and that Eq. ( 1 1 ) can be written in simpler form as 

q(t) = -b(P)T)(t)-a(P)P 

c a (0 = c a (°°; P) + y c ap (°°; P)p(t) + n(t) 

after substituting rjyt) for the deficiency function (and dropping the subscript fi). A broad class of nonlinear 
unsteady responses can be modeled with this model structure where the static term, steady-flow damping term, and 
indicial term parameters can vary nonlinearly with p. The unsteady term has the form of a linear parameter varying 
(LPV) system 2 ' where the parameter is the state, p. instead of a general parameter. In this step, it is assumed that 
each of the four unknown parameters ( C f (oo; /?), C (p (oo; /?), a(/?), b(P) ) can be represented by polynomials in /?. 

The identification problem for the model in Eq. (25) is solved in a two stage iterative process, shown in Fig. 1. 
First model structure determination of the unsteady term polynomials, a(P) and b(P), is performed using stepwise 
regression 2 (SR). To setup the regression equation, a variable y(t) is formed by subtracting the static term and 
approximated steady-flow damping term from the force or moment measurement, CJt), as 


y(t) = C a (0 -C a (co;P)-l C ap (oo; p)p{i ) 

Then the regression equation, showing the discrete measured values at time t(i), is given as 


(26) 


y E (0 = -WPe (0 )1 e (0 + a (Ps (0 )Pe ( 0 ] + £ tl (0 


(27) 


where index E indicates the measured values and £ n (i) is an equation error. Derivatives of measurement, y E , are 
formed using a locally smoothed numerical derivative 24,25 . Parameter estimates are then updated with Maximum 
likelihood (ML) estimation. State and measurement equations for ML estimation are 


id) = ~b(P)n(t) ~ a(P)P E ; >i(t = 0) = ri( 0) 

n E (0 = C aE (0 - C a (oo; p E ( i )) - l - C y (oo; P E ( i))p E (;) + 1/(0, / = 1, 2, ..., N 


(28) 


A second stage identification is done since C Up (oo; P) is only known approximately and the estimates of a(P) 
and b(P) are based on this a priori information. To update the damping term, a second variable is formed 


z(t) = c a (t) - c a (crj; P) - /;( t) 
= y c a p (^;P)p(t) 


(29) 
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If the measured values and new estimates for the unsteady term are substituted into Eq. (29), a new regression 
equation is obtained as 


z /?(0 = y c a p /?( 0 + i = \2,...,n (30) 

The model structure and parameters for the steady-flow damping term can be estimated using Eq. (30). This 
completes the first iteration of the two stage identification. In the next iteration, a(J3) and b((3) can be estimated again 
using the new values and model structure for the steady-flow damping term. This process is continued until 
parameter estimates have converged. 

Step 3. In the last series of dynamic tests, a nonlinear model methodology is required. The last assumption to be 
removed is that superposition of orthogonal planar motions is valid. For this case experiments must be designed to 
produce coupling motions that can be used to test the superposition assumption. If this assumption is proved false 
then experiments designed to allow simultaneous estimation of both longitudinal and lateral-directional general 
model terms are required. From Ref. 3 a set of four experiments are suggested for the case where linear coupling 
occurs and under conditions of nonlinear coupling more complex experiments requiring oscillatory coning may be 
required. These types of experiments were not available for this study. 

V, Experiments 

Dynamic tests were conducted with a 2.5% scale model of the F-16XL (Fig. 2) in the water tunnel at Rolling 
Hills Research Company to identify nonlinear aerodynamic model parameters. For these tests the scale model was 
mounted on a dynamic test rig through a five-component strain-gauge balance (axial force was not measured). The 
dynamic test rig is a computer-controlled system with a sting-mounted double C-strut support system 26 . The 
mounting arrangement rotated the model about the reference center of gravity location of 0.558 c . The tests were 
conducted at a dynamic pressure of 0.81 psf resulting in the flow velocity of 11 inches/sec and a Reynolds number 
of 52x1 0 3 based on the mean aerodynamic chord. Reynolds number values different from that found in flight or 
wind tunnels are acceptable for this study since the data are used to both estimate and predict responses under the 
same conditions. With this relaxed requirement, hydrodynamic flow at the inlet was improved by adding an inlet 
fairing to block flow into the inlet. Flow visualization and measurement equipment inside the model did not allow 
enough space for smooth flow through the body. All data were obtained with the leading edge flaps and trading- 
edge surfaces at 0° deflection. 

Although forced oscillation tests were performed about all three body axes, only roll axis oscillation data are 
considered in this paper. Single-frequency and wide-band (Schroeder sweeps 1517 ) forced-oscillation tests were 
performed about 14 mean angles of attack, a 0 = [10, 15, 20, 25, 30, 32.5, 35, 37.5, 40, 42.5, 45, 47.5, 50, 55, 60], 
with 4 mean offsets in sideslip, j3 n = [0 -5 -10 -15 ], while using 5° amplitude of oscillation. Large amplitude roll 
cases were done at the same a n , for a reduced set of offsets, J3 0 = [0 -5], and amplitudes of [10 15 20 30] degrees. 
Five non-dimensional frequencies, k = [0.066, 0.131, 0.197, 0.262, 0.328], were used for small amplitude tests and 
six non-dimensional frequencies, k = [0.066, 0.095, 0.131, 0.160, 0.197, 0.262], were used for large amplitude roll 
oscillation tests. Data were sampled at 10 Hz with a low pass analog fdter at 5 Hz. For these experiments 30 cycles 
of data were recorded for each inn and then an averaged cycle was formed to minimize measurement noise. The 
averaged measurements were used for small amplitude analysis and 3 cycles of measured data were used for the 
large amplitude analysis. Static data were obtained in 1° increments for angles of attack from 0° to 70° and for 
sideslip angle from -20° to 20°. 


VI. Results and Discussion 

Figure 3a presents the static data for Q (oo;«,/?) as a surface plot and Fig. 3b shows the corresponding contour 
plot. This plot is representative of all the aerodynamic static data in that the surface is relatively smooth and linear at 
low and very high angles of attack. In the mid angle-of-attack range, 20° < a< 50°, fairly severe static nonlinearities 
occur and the model is slightly asymmetric. This may limit the effectiveness of a simple polynomial representation 
of the underlying static nonlinearities as part of the overall model. 
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Harmonic analysis of the small-amplitude oscillatory data was performed on single oscillations formed by 
averaging 30 cycles of data. This analysis was done for each of the 5 frequencies. Figure 4 shows representative 
input / output data used in this analysis. The figure shows averaged small-amplitude roll oscillation time histories for 
P(t), p(t), and C ( (t;a(t),/3(t)) at a low and high frequency with offset J3 0 = 0° and a 0 = 30°. With small-amplitude 
oscillations a linear model is generally expected to be adequate for all a 0 and p n set points. Harmonic analysis 
confirmed this for most conditions tested, however, at a 0 = 42.5°, Po = 0°, k = 0.066, and at a 0 = [37.5°, 40°], with 
Po = -5°, and k = 0.066, the Multiple Correlation Coefficient, R 2 was relatively low (less than 0.6), indicating less 
adequate linear models for these conditions. 

Analysis of the small amplitude data provides a basis for the nonlinear analysis. The variation of linear model 
terms with angle of attack and sideslip defines the initial model structure used in the nonlinear analysis. First-order 
harmonic analysis provided in-phase and out-of-phase components, ( C f p , C fp ). Out-of-phase ( C fp ), damping terms 

are shown in Fig. 5 over the range of angle of attack, frequencies, and p 0 offsets. These results show a strong 
frequency dependence or unsteady behavior for 25° < a < 40°, with Po = 0°, and for 25° < a < 37.5°, with p 0 = -5°. 
For offsets of p 0 = [-10°, -15°] the unsteady behavior is practically nonexistent. Overall a linear model is reasonable 
for the small amplitude data, thus Two-Step Regression 10 ' 27 was applied to the in-phase and out-of-phase data to 
produce a linear unsteady model where appropriate. As an example, Fig. 6 highlights the linear model parameters 
( Q/?(°°) , C fp (co) , a, rj) over the range of angle of attack for p Q = 0. At a 0 = 42.5°, the parameter estimates are in 

doubt due to poor R 2 value. For this example, C^(oo), estimated by regression matches well with the slopes 

estimated directly from the static data indicating relatively good identification. The steady-flow damping term, 
Cf p (oo) , is relatively near -0.2, except at a 0 = [35°, 37.5°] where it moves closer to 0. The unsteady term gain, a/b, 

confirms that for 25° < a < 40°, unsteady behavior is present. The gain values become very small outside of this 
range of angle of attack indicating no significant unsteady behavior. The non-dimensional time constant, zy, has a 
maximum value of 20 in the unsteady region and values less 
than 10 outside this region. The lack of unsteady behavior at 
do = 50° resulted in a poor zy estimate with large error 
bounds. Figure 7 summarizes C (p { oo) for each /?„ offset at a n 

= [30°, 37.5°]. For a 0 = 37.5°, p 0 = -5°, Q p (oo) is rejected 

due to poor R 2 and the lack of consistent frequency 
distribution in the harmonic coefficients. These estimates 
have large error bounds reflecting the lack of information 
content. The lack of unsteady behavior and inconsistent 
frequency distribution in this range can be seen in Fig. 5 and 
Fig. 8. A restriction for the linear model in Eq. (14) is that it 
only provides an adequate model for portions of the envelope 
where unsteady behavior occurs, as indicated by region II in 
Fig. 8. 

Figures 9a and 9b show sample time histories of the large 
amplitude oscillatory data at a 0 = 30° and 37.5°, each case 
with Po = 0°. The figures show four cycles of oscillation for 
each of the four frequency cases used in analysis, k = [0.066, 

0.095, 0.131, 0.160]. The time histories have been 
concatenated to facilitate data handling. The first cycle of predicted response is eliminated to remove any initial 
transients and ensure that steady harmonic motion is used in calculating residuals. The amplitude, (f>A, is the same for 
each case shown, however, the resulting a and P oscillations change during roll oscillations according to Eq. (24). 
The distortion of the rolling moment response from a steady harmonic wave, in particular for the low frequency 
cases, indicates nonlinear behavior. 

Harmonic analysis was applied to the large amplitude data. First in-phase and out-of-phase components were 
calculated from the 1 st harmonic coefficients. Plots of the out-of-phase components over the range of angle of attack 
and frequencies are shown in Fig. 10 with Po = 0. This graph shows the effect of oscillation amplitude on the 
damping coefficient relative to the small amplitude case shown in Fig. 5. For the large amplitude case the unsteady 
behavior is observed over approximately the same range of angle of attack (25° < a < 40°) as the low amplitude 
case, however, the magnitude of the out-of-phase damping is reduced by almost half. Figure 11 shows the multiple 


Table I. Linear model parameter estimates 
and standard errors (in parentheses) for C f , 
using two identification methods, at a 0 = 30°, 
Po = 0°, fa = 30°. 


Parameter 

2-Step 

Analysis 

Relaxation 

Method 

Q/j(°°) 

0.06 

(0.012) 

0.063 

(n/a) 

Ctp(P°) 

-0.28 

(0.010) 

-0.321 

(0.0026) 

a 

0.305 

(0.0049) 

0.309 

(0.001) 

b 

0.15 

(0.0025) 

0.156 

(0.0009) 

Tl 

14.2 

(0.24) 

14.5 

(0.0020) 

R 2 

0.99 

0.98 
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correlation coefficient, R 2 , from the large amplitude analysis. The upper graphic shows the variation of R 2 with 
angle of attack and non-dimensional frequency, k. In this case the low frequency run with k=0.066, has the lowest 
values of R 2 for 35° < a< 40°, indicating some degradation in the linear model at these angles of attack. The lower 
graphic shows the variation of R 2 with angle of attack and harmonic model order for this worst case frequency, 
k=0.066. From the harmonic analysis, only 
adding a 3 rd order harmonic improved the 
result. One explanation of this could be the 
presence of a 3 rd order polynomial nonlinearity 
in either the static or damping term. In order to 
identify the model in this case the relaxation 
method is applied. 

To demonstrate the method, two examples 
are considered using the data in Fig. 9. The first 
example, at a 0 = 30°, considers linear unsteady 
behavior that can be adequately modeled by the 
linear model methodology. Application of the 
relaxation method allows for more complex 
modeling but in this case it is used to confirm 
the linear analysis and validate the method. 

Parameter estimates and their standard errors 
using both methods are given in Table I. Both 
methods provide nominally the same models. 

Figure 12 shows a comparison of measured and predicted response, in a phase-plane plot, using the linear model 
obtained from the relaxation method. Further refinements are possible, such as modeling the static data, shown in 
Fig. 13, with a cubic or even 5 th order polynomial to more exactly capture the true nature of the system. However, in 
this case the static moment contribution is a factor of 10 smaller than the damping term contribution; consequently 
the increased model complexity is not justified. R 2 values in Fig. 1 1 confirm this assumption for all frequencies. 

Measurements at = 37.5° provide an example of response with nonlinear behavior, based on R 2 values shown 
in figure 11, and strong unsteady behavior based on frequency dependence of in-phase and out-of-phase coefficients 
shown in figure 10. Figure 9b shows the responses at all frequencies are distorted from a linear sinusoidal response. 
Since the harmonic analysis inferred that at least a cubic nonlinearity is present and the static component is a 
relatively large component in this case, a cubic 
representation of the static data was initially 
implemented. It was found, however, that a 5 th 
order polynomial was required for an adequate 
static model. Static rolling-moment coefficient, 

Q(oo ;« 0 = 37.5,/?) , is presented in Fig. 14 

together with its approximation by the 5 th order 
polynomial. Results from linear regression 
provided initial estimates for the damping and 
unsteady terms in the postulated nonlinear 
model. Table II shows linear models estimated 
from both small and large amplitude data. 

Linear analysis of the large amplitude data 
produced a slightly different model compared 
to the small amplitude case reflecting the 
amplitude dependence of nonlinear systems. 

Estimated parameters are shown with their 
standard errors below the estimate in parentheses. Table III shows the final estimated model from the relaxation 
method applied at a n = 37.5° and the corresponding overall R 2 for all four frequencies. The parameter covariance 
matrix for this case contained some large pair-wise correlations indicating identifiability may be a problem. Figure 
15 shows the measured and predicted rolling moments for the four frequencies considered. These graphs all show 
distortion from the regular elliptic shape associated with linear response. For each frequency, the nonlinear model 
prediction matches the response well. At the lowest frequency, k = 0.066, the distortion or nonlinearity is greatest, 
but the model match is still very good. Although the fit is very good the model is not completely satisfactory due to 
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Table III. Nonlinear model for Q at a 0 = 37.5°, p 0 = 0°. 


Parameter 

Nonlinear Model 

QC °o;y?) 

0.007 +0.317,0 +0.202/-12.786/ 
-2.371/ +91.123/ 

C, p (^,P) 

-0.786 -1.791/+97.616/ 
(0.0045) (0.34) (4.48) 

a(/3) 

0.527 -43.668P 2 -507.228/ 
(0.0051) (0.39) (5.096) 

b(P) 

0.624 -8.109/ 
(0.011) (0.14) 

II 

o 

3.63 

(0.027) 

R 2 

0.97 


Table II. Linear models for C f at a 0 = 37.5°, p 0 = 0°. 


Parameter 

Linear Model 

Linear Model 


<t>A = 5° 

<t> A - 30° 

c v (°°) 

0.15 

-0.015 


(0.045) 

(0.0096) 

Q/co;/?) 

-0.13 

-0.31 


(0.049) 

(0.020) 

a(P) 

0.47 

0.242 


(0.023) 

(0.0095) 

b(P) 

0.28 

0.44 


(0.032) 

(0.035) 

hiP) 

8.0 

5.1 


(0.90) 

(0.41) 

R 2 

0.90 (k=0.066) 

0.78 (k=0.066) 





the inconsistency with the small-amplitude linear analysis shown in Table II. If the nonlinear model can represent 
the linear case as a subset then, for example, at f} = 0, the nonlinear model estimate, C fp (oo; /? = ()), should be close 

to the corresponding linear model analysis of the small amplitude data, a value of -0.13. The differences may reflect 
a number of possible sources of error besides the model structure for the damping and unsteady terms, such as data 
accuracy and information content of the data. 

VII. Concluding Remarks 

This paper introduces a general formulation of the aerodynamic model for aircraft that includes nonlinear 
unsteady aerodynamics and suggests a progressive series of steps toward development of a nonlinear model 
identification methodology that defines the mathematical model, algorithms for identification, and required test 
techniques. The mathematical model presented can predict a broad class of aerodynamic responses of rigid-body 
aircraft. An important element of the methodology is the form used for the aerodynamic model. This form retains 
conventional static and rotary aerodynamic terms that have traditionally provided substantial engineering 
information to the flight dynamics community. This structure allows easy interpretation of the model parameters by 
retaining conventional stability and control derivatives for static and dynamic terms. Unsteady terms in the form of 
indicial functions are added to expand the model structure and allow frequency dependent phenomena to be 
modeled. In the nonlinear case, parameters in each of these terms may vary with angle of attack and sideslip. Thus, 
static and rotaiy dynamic terms can be represented as polynomials or splines and the indicial term as a linear 
parameter varying system. This approach produces a straightforward and efficient model structure for simulation 
and design since these models are more general than conventional locally-linearized models but still easily 
incorporated into conventional table look up formats. 

As a practical matter, it is desirable to develop an adequate model that represents the physics of interest without 
unnecessary complexity. In addition, the mathematical model must lend itself to identification from experimental 
data. This suggests model development should start with the least complicated model and add complexity only as 
needed. Initial steps for developing a nonlinear model identification methodology with this philosophy are applied to 
static and forced oscillation data to identify rolling moment coefficient. For nonlinear model identification a two- 
stage identification procedure was used. From experimental data it was assumed that the static and rotary-dynamic 
terms were known. The method combined stepwise regression and maximum likelihood estimation to determine 
model structure and parameter estimates for the indicial term. When a priori information about the rotary term is 
poor, a two-stage optimization algorithm can be used to identify both the unsteady and rotaiy terms. The method 
was applied to a linear unsteady case that validated the approach against a 2-Step Regression method. A second case 
considered a nonlinear unsteady model that provided good fit to the nonlinear dynamic response measurements. The 
nonlinear model was not completely satisfactory, however, because it did not reduce to the linear model for small- 
amplitude motions. High correlations for some parameter estimates suggest an identifiability issue that requires 
further review of the methodology, experiment design, and measurement accuracy. Nonlinear modeling 
experiments, in particular, require careful attention to these issues. An assumption in the proposed candidate 
methodology that new measurements could be formed by differencing the measured rolling moment with the static 
and damping terms to form a new indicial measurement was not validated by this study. This approach assumes that 
the new measurement will be corrupted primarily by measurement noise and not model error. The successful 
application of this assumption depended on the relative magnitudes of the three terms in the model. 

Further development of nonlinear model identification methodology is needed to realize the capability for 
predicting aircraft aerodynamic response under arbitrary motion. Future work should refine experiment design and 
identification algorithms, and investigate a range of input amplitudes and wide-band inputs for maximizing 
information content and testing efficiency. Advanced test techniques such as oscillatory coning, and planar forced- 
oscillation in heave and sway should be studied since these offer a direct measurement of unsteady acceleration 
terms and can provide validation of the nonlinear model identification methodology. 
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two stage optimization 



Eqn. A: z E (t) = \Jfj C a p p(0 

Eqn-B: y E (t) = C fl (f)-C fl (co;a 0 ,£)-^jc (oo;a 0 ,^) p(#) 

Figure 1. Block diagram of model identification procedure using stepwise regression 
(SR) and a maximum likelihood (ML) estimation. 



Figure 2. Three-view drawing of 2.5% F-16XL model. 


Figure 3a. Surface plot of C t (oo; a, J3) . 
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Figure 3b. Contour plot of Q(co; a,/3). 



a, deg 


Figure 5. Variation of out-of-phase component with 
angle of attack, offsets in sideslip, and frequency of roll 
oscillations. 
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Figure 7. Variation of damping-in-roll parameter with 
sideslip angle at two values of angle of attack. 


Figure 6. Estimated linear model parameters for rolling- 
moment coefficient using 2-Step Regression. 



a 0 , deg 


Figure 8. Three regions with different dynamics 
for roll moment in-phase and out-of-phase 
components. 
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Figure 9a. Time histories of the large amplitude oscillatory data at a 0 = 30°, p 0 
0°, and four frequencies, k = [0.066, 0.095, 0.131, 0.160]. 



time, s 

Figure 9b. Time histories of the large amplitude oscillatory data at a 0 = 
37.5°, p 0 = 0°, and four frequencies, k = [0.066, 0.095, 0.131, 0.160]. 
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Figure 10. Variation of first-order harmonic 
coefficient with angle of attack and frequency. 



Figure 11. Variation of Multiple Corelation Coefficient 
with angle of attack, frequency, and order of harmonic 
coefficient. 
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Figure 13. Comparison of measured and 
computed rolling-moment coefficient from 
3 rd order polynomical at a 0 = 30°. 



Figure 12. Comparison of measured and computed rolling- 
moment coefficient from relaxation method at a 0 = 30°. 



p, deg 

Figure 14. Comparison of measured and 
computed rolling-moment coefficient from 
5 th order polynomial at a 0 = 37.5°. 





Figure 15. Comparison of measured and computed rolling- 
moment coefficient from relaxation method at a 0 = 37.5°. 
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